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We propose an initialization procedure for the density-matrix renormalization group 
(DMRG): the recursive sweep method. In a conventional DMRG calculation, the infinite- 
algorithm, where two new sites are added to the system at each step, has been used to reach 
the target system size. We then need to obtain the ground state for a different system size 
for every site addition, so 1) it is difficult to supply a good initial vector for the numerical 
diagonalization for the ground state, and 2) when the system reduced to a ID system consists 
of an array of nonequivalent sites as in ladders or Hubbard-Holstein model, special care has 
to be taken. Our procedure, which we call the recursive sweep method, provides a solution to 
these problems and in fact provides a faster algorithm for the Hubbard model as well as more 
complicated ones such as the Hubbard-Holstein model. 

KEYWORDS: DMRG, infinite algorithm DMRG, finite algorithm DMRG, exact diagonalization, ladder, 
Hubbard model, electron-phonon interaction, phonon, Hubbard-Holstein model 



Introduction After fifteen years since the density- 
matrix renormalization group (DMRG) ' was intro- 
duced as a new numerical method to obtain ground- 
and low-energy excited states for one-dimensional many- 
body systems, DMRG continues to enjoy its status as 
one of the most powerful numerical tools in condensed 
matter physics. It has been successfully applied to study 
various problems, 3 which include correlated electron sys- 
tems such as the Hubbard model 4-6 and the t-J model, 7 
quantum systems at finite temperatures, 8-10 quantum 
chemistry, 11 quantum Hall systems, 12 dynamical quan- 
tities. 13-15 The extreme accuracy obtained by DMRG 
in many cases has attracted interests for its mathemat- 
ical reasons, and such viewpoint is also pursued vigor- 
ously (see, e.g., 16 ' 17 ). One of the recent breakthroughs 
in the field is the development of several new methods to 
apply DMRG to real-time simulations of quantum sys- 
tems. 18-20 Several review articles 21-23 describe the his- 
tory of DMRG and its wide range of application. 

In DMRG, we consider a system consisting of sites 
aligned in a one-dimensional open-boundary chain. We 
systematically transform the Hilbert space of the sys- 
tem to one whose dimension is drastically reduced from 
the original one, while keeping the low-energy physics as 
unchanged as possible. In order to do this, we separate 
the system into parts, which we call blocks or subchains, 
and reduce the number of basis functions to describe 
the blocks, which we call the basis size for the blocks. 
What we actually do is to choose a basis that is fit to 
express one (or a few) quantum-mechanical state(s) of 
the original system (target state(s)), and expect that the 
low-energy excitations are recaptured in the chosen ba- 
sis. The basis size to express physical states is reduced 
with the use of a partial density matrix calculated for 
some parts of a superblock, which is the whole or a par- 
tial chain and is usually constructed from two blocks and 
two sites. The part to calculate the partial density ma- 
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Fig. 1. DMRG calculation. The one-dimensional lattice model is 
treated by iteratively adding sites to two subchains to construct 
supcrblocks from two subchains and two sites between them. In 
enlarging a subchain to contain one more site, the partial den- 
sity matrix for the new block is calculated for the target state of 
the superblock. The eigenstates with m largest eigenvalues are 
retained as the new basis for the enlarged subchain. The trunca- 
tion error is assessed by calculating the sum r\ of the eigenvalues 
for the discarded eigenstates. In the new basis, the matrix rep- 
resentations of the Hamiltonian and operators for the new block 
are calculated and stored for later use. 



trix is separated at some cut point from the rest of the 
system, which is called the environment block. While the 
number of states to be considered in an exact diagonal- 
ization increases exponentially with the number of sites, 
the contribution of high-energy configurations to the low- 
energy physics is usually exponentially small. With the 
use of the DMRG algorithm, most important states can 
be systematically found. In many cases, the dimension 
of the Hilbert space can be kept much smaller than in 
the exact diagonalization, and almost exact results can 
be obtained with moderate computational efforts. 

In DMRG we need to start with a warm-up process, 
which is the process to iteratively add new sites to one 
or more subchains, until the Hilbert space of the sys- 
tem with the desired number of sites can be represented 
with the reduced basis. If the system consists of homoge- 
neous sites, we can just add two new sites between two 
subchains of the same length I to make a superblock of 
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21 + 2 sites in order to obtain a subchain that has I + 1 
sites, and repeat this until the superblock acquires the 
desired number of sites. This process is called the infinite 
algorithm DMRG; we can repeat the addition of a site to 
the subchain indefinitely. Once this warm-up procedure 
has finished, we can enhance the quality of the basis by 
systematically moving the cut location back and forth 
within the chain (the finite algorithm DMRG). 

In this letter, we focus on the warm-up process. When 
we increase the number of sites by two in the infinite- 
algorithm, it is in general difficult to supply a good initial 
vector for the numerical diagonalization for the target 
state, so usually a random vector is plugged. Our new 
procedure reduces the number of numerical diagonaliza- 
tion from a random vector, with the use of the finite 
algorithm DMRG. 

Method We illustrate the new procedure for the case 
of N c original sites per full site. Our idea here is to modify 
the infinite algorithm in the conventional DMRG, which 
adds two sites at the center of the superblock at each 
step, so that two full sites are added at the center of 
the superblock per cycle. As we elaborate in the below, 
for each of the left and right subchains, we retain N c 
subchains that have respectively 0, 1, . . . , N c — 1 original 
sites besides n full sites, instead of just one subchain 
each. The longest subchain among the left subchains and 
the longest among the right ones are used to construct a 
new superblock, and with the use of the finite algorithm, 
those subchains are enlarged by 1, 2, . . . , N c original sites. 
Then we can repeat the same process for a superblock 
that is two full sites longer. 

Now we explain how this can be done. We align the 
original sites in a one-dimensional chain so that the same 
kind of sites appear periodically along the chain, as 

• • • ) — (b b^-i) — (b ^n.-i) — (•••, 

where the original sites are denoted as bo, . . . ,bjv _i, and 
each full site is enclosed in (•••). 
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Fig. 2. The recursive sweep method illustrated for the case of 
n = 2 and N c = 3. A full site, s, is represented as an oval, and an 
original site, b^, as a circle. An arrow shows how an original site 
is added to a subchain, which is represented as a round cornered 
rectangle. The first steps to enlarge the left and right subchains 
are marked as (i) and (ii), as in the text. For the definition of 
the superblocks Ag°' ±1,±2 '' and Ag°\ see the text. 



We denote the full site as s, and a sequence of n full 
sites as s". Suppose we have the left and right subchains 



L nt i and R n ,i {i = 0, 1, . . . , N c — 1) respectively, 

L n ,o = s", R n ,o = s", 

L n ,i — s" — bo, R n ,i — bjv c _i — s™, 

L n ,2 = s" — bo — bi, R n ^ — bjy c -2 — bjv c -i — s", 

L n ,N e -i = s" — b bjv _i, and 

Rn,N c -l — bi bjy c _2 — bjv c _i — S™, 

where L n> i (R n ,i) has (a) n full sites from the left(right) 
edge of the original chain and (b) i(— 0, 1, . . . , N c — 1) 
original sites from the (n+I)-th site s n (sl-u-i)- We en- 
large L n! jv c _i and R n ^ a -i by N c original sites to obtain 
2N C subchains L n+ i_, and R n +\^ (i — 0, 1, ... , N c — I) 
as follows. 

First, we construct a superblock A 2 °^ +2 that has 2ra + 2 
full sites: 

is constructed from L„jv c _i, two sites p a and pp, and 
Rn,N c -i, aligned in this order from the left. We calcu- 
late the target state(s) for A 2 ^ +2 . Then we can obtain 
(i) L n+ ifi by calculating the partial density matrix for 
L n ,N c -i and p a combined, and in the same way, (ii) 
i?„+i ; o from and R n ^ a -i combined. We write this 
as 

(i) [L n ,N e -x — bjv c _i(p a )] — » L n+ i :0 , 

and 

(n)Rn +1 fl <— [b (p/3) — R n ,N c -i]- 

Next, as in the finite algorithm, we move the cut lo- 
cation to the right by one original site. This is possible 
because we have i?„ i jv t ,_2. The new superblock, A 2 | l '' +2 
also has 2n + 2 full sites: 

A 2n+2 = ^n+1,0— bofoo) — bi(p 7 )— Rn,N c -2 

is constructed from L„ + i.o, two original sites pp and p 7 , 
and Rn.Nc-2- Then we obtain L„ + i.i by calculating the 
partial density matrix for L n+ i^ and pp combined: 

[£n+i,o — b (p/3)] — > 

This can be repeated N c — 2 more times, where for each i 
(i = 2, . . . , iV c — 1), we obtain L n +ij from the superblock 

A-2n+2 that is constructed from L n+ i two original 
sites, and R n ,N c -i-i- 

^■2n+2 = bj_l bi RnM a -l-i, 

[L n +l,i-l bj_i] — » L n+ i^. 

Also, by moving the cut location to the left N c — 1 
times, we make superblocks A 2 ^ 2 (i = 1, . . . , N c — 1) 
that consists of L„ + i jv c -i-ij two original sites and 
Rn,i-u to obtain R n +i,t (i = 1, ■ • • , N c — 1): 

^2n+2 ~ ^n,iV c -1-t bN c -l-i — bjv e -t •Rn+l.i-lj 

Rn+1A <— [b Nc _ i Rn+l,i-l], 

for % = l,...,N c - 1. 

In these steps, the initial vector for each of the exact 
diagonalizations can be calculated from the target state 
obtained in the most recent step (for A 2n ^ 2 , we use the 
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target state for A^!^, rather than A^+2 )• The cal- 
culations for {L n+ \^i\ and can be done without 
data dependence to each other in a parallel computer. 

Now we have L n+ i^ (i = 0, . . . , N c — 1) and R n +i,i (i = 
0, . . . , iV c — 1), so we have enlarged the subchains by one 
full site. We have to obtain the target states from random 
vectors only once. This procedure can be repeated for the 
incremented values of n, until the superblock has all the 
L full sites. 

Note that, because DMRG limits the number of the 
basis functions to to to describe finite-size blocks, a cal- 
culation with any choice of basis should converge to the 
result of the exact diagonalization in the to — > oo limit. 
The questions here are: 

• how fast the calculation converges, and 

• whether we can calculate with enough accuracy 
within a practical computation time. 

We also note that the recursive sweep method de- 
scribed here is different from (a) the initialization by 
Liang and Pang 24 for a rectangle lattice where the ra- 
tio of horizontal and vertical hopping parameters is con- 
trolled through several finite-algorithm sweeps that in- 
volve the whole system, (b) the warm-up with a conven- 
tional numerical renormalization group 25 by Xiang, 5 the 
two-step DMRG by Moukouri and Caron, 26 and (c) the 
quantum information entropy-based approach by Legeza 
and Solyom. 27 

Choice of targets Because environment blocks are 
shorter than the enlarged block in the finite-algorithm 
steps, it is advisable to target several states having dif- 
ferent quantum numbers {e.g., numbers of up-spin elec- 
trons and down-spin electrons) near the target state. 
In calculating for the half-filled, L-site system (L: an 
even integer), we target not only the ground state for 
(n,f,ni) = (L/2,L/2) but also the ground states for 
(n T + d h] ,n L +d Lr ) {dj,d f = 0,±1). 
Results 

The Hubbard model We first apply the recursive 
sweep method to the half-filled Hubbard chain with ho- 
mogeneous sites, 

H = -J2 M^+i.a + H.c.) + ]T U(n lA - l)(n a - 1), 

2,(7 1 

(1) 

by considering multiple sites as one composite full site. 
In each iteration of the recursive sweep method, we add 
two composite sites to the superblock. When the length 
of the superblock equals or exceeds the desired length 
Lq, we can stop the recursive sweep method and start 
the finite-algorithm sweeps. 

As shown in Figure 3, the time required for the warm- 
up stage is almost halved in the case of the 84-site Hub- 
bard model with U/t = 2 by the use of this method. The 
longer cycle we take, generally the faster the calculation 
becomes. We also plot the ground state energy per site 
in Figure 3. The ground state energy becomes slightly 
worse than in the conventional infinite-algorithm warm- 
up with the same number of retained states to = 50. 
When we increase m to m = 80, however, the energy be- 
comes better than in the infinite-algorithm with to = 50 



03 

E 

Q. 

O 



£ 
>. 

i— 

CD 
C 

LU 



1000 
750 
500 
250 


-1.3365 
-1.337 

-1.3375 
-1.338 



m=80 



m=50 

conventional method 



—El 




m=800 (Finite algorithm) 



1 2 4 6 8 10 12 14 
Number of sites / full site 

Fig. 3. Upper panel: The CPU time for the warm-up for the 84- 
site Hubbard model with U/t = 2, for the number of retained 
eigcnstates m = 30, 50, 80, as function of the number of sites 
in a full site for the recursive sweep method. The calculations 
were conducted on a workstation with two Xeon 5160 chips and 
4GB of RAM. Lower panel: The ground-state energy per site 
at the end of the warm-up for the 84-site Hubbard model with 
U/t = 2, for the number of retained eigenstates m = 50, 80, as 
function of the number of sites in a full site for the recursive 
sweep method. The result for the conventional infinite method 
is also shown. The true energy, which can be estimated with the 
finite-algorithm DMRG with a large m = 800, is -1.337 785 42 
and shown as a dotted line. Other lines are guide to the eye. 
The ground state energy per site for the infinite-length Hubbard 
model with U = 2, obtained with the exact solution by Lieb and 
Wu, is -1/2 - J (x)J 1 (x)/x(l + eW)dx = -1.344 374 34. 



for up to seven sites per cycle, while the CPU time spent 
is still significantly lower. So a small increase in to can 
compensate the increased error due to the use of the re- 
cursive sweep algorithm with N c = 2 — 7. 

The Hubbard- Holstein model When the system is a 
repetition of three or more types of sites, as in the case 
of the pseudo-site method for the Holstein phonons 28 or 
multi-leg ladder systems, the choice of the warm-up pro- 
cess in DMRG is not obvious. Hereafter, we call one orig- 
inal site, which becomes the period of the repetition of 
the (pseudo-)sites, as a full site, (a) If we enlarge the 
two subchains in a mirror-symmetric way, at most of the 
steps, we have two incomplete full sites at the center of 
the superblock. Then the environment for the new sites 
at the center is much different from the later stages of 
the calculation, when the full site is completely within 
the superblock. (b) Alternatively, we can increase the 
length of the subchain that starts from one edge of the 
chain by iteratively using short subchains whose whole 
basis can be exactly treated, in a way in which the sites 
in the superblock constitute a number of full sites at each 
step. In this case the new sites are always added near the 
edge of the chain, even when they are around the center 
of the original system. So both of the approaches (a)(b) 
have shortcomings that should deteriorate the choice of 
the reduced basis, which can be overcome with the use 
of the recursive sweep method as follows. 

Here we consider the Hubbard-Holstein model 29 on a 
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one-dimensional chain, 

i,(7 i 



(2) 



Here, Ci tCy annihilates an electron with spin er(=j,J,) at 
site i, ni t „ = c\ a ci. a is the electron number, g is the 
clcctron-phonon coupling, and is the phonon anni- 
hilator at site i. For the application of the pseudo-site 
method to the Hubbard-Holstein model, the present au- 
thor, in a collaboration with Arita and Aoki, 30 has pre- 
viously come up with another method, the compensa- 
tion method, to make improve the choice of basis for the 
infinite-algorithm warm-up. Because DMRG is a vari- 
ational method, the calculated ground state energy is 
always higher than the actual value; the lower energy 
means the better convergence. 

As we show in Figure 4, the recursive sweep method 
warm-up (solid circle in the plot) results in a much better 
energy per site that does not fluctuate and continues to 
decrease as the number of full sites in the left subchain 
is increased. Similar results have also been observed for 
other parameter sets. 
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Fig. 4. The energy per electron in units of t against the length 
of the chain at intermediate stages of the warm-up process, cal- 
culated (i) with the infinite algorithm without the compensa- 
tion method, (ii) with the infinite algorithm with the compen- 
sation method, and (iii) with the recursive sweep method, for 
the Hubbard-Holstein model with (U/t, g/t,LOo/t) = (2,3,5). 
2^=16 is chosen for the cutoff in the phonon number. The size 
of the basis is m = 70 for all the calculations. The memory 
required for the calculation is around 700MB for the recursive 
sweep method, when all subchains and operators for calculation 
of various correlation functions are retained for later use. 



Summary The recursive sweep algorithm allows an 
initialization of a DMRG calculation when the system is 
a repetition of multiple kinds of sites, without the need of 



treating incomplete cycles of sites or the need of adding 
new sites at locations much off the center of a superblock. 
It is an extention of the infinite algorithm DMRG by 
the finite algorithm sweeps, over one cycle of added sites 
to both directions along the chain, whose computational 
cost is much smaller than in the conventional infinite 
algorithm. We have demonstrated that the present al- 
gorithm improves the calculation for a system of local 
phonons coupled to correlated electrons compared to our 
previous method. The recursive sweep algorithm is also 
applicable to a chain with homogeneous sites. We have 
demonstrated that we can both reduce the calculation 
time and improve the energy when we apply the algo- 
rithm with a small increase in the number of retained 
states. 
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